rm(list=ls())

load("TSCS_data.RData")

# Table E.4

lm_down_1.1 <- plm(downturn ~ lag(downturn,1:2)  + 
                     antipluralism_gov_seat_share_lag + 
                     citizen_support_lag2 +
                     lag(polarization,1) +
                     lag(clientelism_gov_seat_share,1) +
                     lag(region_citizen_support,1) + 
                     lag(region_antipluralism_gov_head_seat_share,1) + 
                     lag(region_liberal_democracy,1) + 
                     lag(democratic_stock,1) + 
                     lag(GDP_capita,1) + 
                     lag(life_expectancy,1),
                   data = panel,
                   model = "within",
                   effect = "twoways")

lm_down_1.2 <- plm(downturn ~ lag(downturn,1:2)  +
                     antipluralism_gov_seat_share_lag*citizen_support_lag2 + 
                     lag(polarization,1) +
                     lag(clientelism_gov_seat_share,1) + 
                     lag(region_citizen_support,1) + 
                     lag(region_antipluralism_gov_head_seat_share,1) + 
                     lag(region_liberal_democracy,1) + 
                     lag(democratic_stock,1) + 
                     lag(GDP_capita,1) + 
                     lag(life_expectancy,1),
                   data = panel,
                   model = "within",
                   effect = "twoways")

cmodel_down_1.1 <-
  pgmm(
    downturn ~ lag(downturn,1:2) + 
      antipluralism_gov_seat_share_lag +
      citizen_support_lag2 +
      lag(polarization,1) + 
      lag(clientelism_gov_seat_share,1) + 
      lag(region_citizen_support,1) + 
      lag(region_antipluralism_gov_head_seat_share,1) + 
      lag(region_liberal_democracy,1) + 
      lag(democratic_stock,1) +
      lag(GDP_capita,1) + 
      lag(life_expectancy,1) |
      lag(downturn, 3:5),
    panel,
    effect = "individual",
    model = "onestep",
    transformation = "ld",
    indexes = c("country_name", "year"),
    robust = T
  )
 
cmodel_down_1.2 <-
  pgmm(
    downturn ~ lag(downturn,1:2)  +
      antipluralism_gov_seat_share_lag*citizen_support_lag2 +
      lag(polarization,1) + 
      lag(clientelism_gov_seat_share,1) + 
      lag(region_citizen_support,1) + 
      lag(region_antipluralism_gov_head_seat_share,1) + 
      lag(region_liberal_democracy,1) + 
      lag(democratic_stock,1) + 
      lag(GDP_capita,1) + 
      lag(life_expectancy,1) |
      lag(downturn, 3:5),
    panel,
    effect = "individual",
    model = "onestep",
    transformation = "ld",
    indexes = c("country_name", "year"),
    robust = T)

summary(cmodel_down_1.2)

# SEs and statistical tests for models

modslin <- list(lm_down_1.1,lm_down_1.2)
vcm_lm <- lapply(modslin, function(x) plm::vcovBK(x, cluster=c("group")))
rse_lm <- lapply(vcm_lm, function(x) sqrt(diag(x)))
pval_lm <- lapply(1:length(modslin), function(i) {
  pt(abs(coef(modslin[[i]]) / rse_lm[[i]]), df = 1608, lower.tail = FALSE) * 2
})
model1_wool <- sapply(modslin, function(x) round(pbgtest(x, order=1)$p.value, 3))[c(1,0)] # Wooldridge test for serial correlation
model2_wool <- sapply(modslin, function(x) round(pbgtest(x, order=1)$p.value, 3))[c(2,0)] # Wooldridge test for serial correlation

modsmom <- list(cmodel_down_1.1,cmodel_down_1.2)
vcm_gmm <- lapply(modsmom, function(x) plm::vcovHC(x, cluster=c("group")))
rse_gmm <- lapply(vcm_gmm, function(x) sqrt(diag(x)))
pval_gmm <- lapply(1:length(modsmom), function(i) {
  pt(abs(coef(modsmom[[i]]) / rse_gmm[[i]]), df = 1608, lower.tail = FALSE) * 2
})
model1_AB <- sapply(modsmom, function(x) round(mtest(x, order=2, vcov=plm::vcovHC(x, cluster="group"))[["p.value"]][[1]] ,3))[c(1,0)] # AR(2) p-val
model2_AB <- sapply(modsmom, function(x) round(mtest(x, order=2, vcov=plm::vcovHC(x, cluster="group"))[["p.value"]][[1]] ,3))[c(2,0)] # AR(2) p-val
model1_H <- sapply(modsmom, function(x) round(sargan(x, weights="twosteps")[["p.value"]][[1]], 3))[c(1,0)] # Hansen test p-value
model2_H <- sapply(modsmom, function(x) round(sargan(x, weights="twosteps")[["p.value"]][[1]], 3))[c(2,0)] # Hansen test p-value

stargazer(modslin,modsmom,
          covariate.labels=c("Backsliding t-1",
                             "Backsliding t-2", 
                             "Anti-pluralist gov. t-1",
                             "Citizen Support t-1",
                             "Anti-pluralist gov. t-1 x Citizen support t-2",
                             "Polarization t-1",
                             "Clientelism t-1",
                             "Region citizen support t-1",
                             "Region anti-pl. government t-1",
                             "Region democracy t-1",
                             "Democratic stock t-1",
                             "log(GDP t-1)",
                             "Life expectancy t-1"),
          order = c(1,2,3,4,13,5,6,7,8,9,10,11,12),
          type="latex",
          se = list(rse_lm[[1]],rse_lm[[2]],rse_gmm[[1]],rse_gmm[[2]]),
          p = list(pval_lm[[1]],pval_lm[[2]],pval_gmm[[1]],pval_gmm[[2]]),
          omit.stat	= c("f","n","adj.rsq"),
          dep.var.caption	= "",
          dep.var.labels	= "Backsliding",
          add.lines=list(c("N Observations","1796","1796","3493","3493"),
                         c("N Countries","99","99","100","100"),
                         c("Wooldridge test (p-value)", model1_wool, model2_wool,"",""),
                         c("Arellano-Bond test (p-value)","","",model1_AB,model2_AB),
                         c("Hansen test (p-value) ", "", "",model1_H,model2_H)
          ),
          no.space = T,
          star.cutoffs = c(0.1,0.05,0.01),
          out = "tab_e4.tex"
)
